.. currentmodule:: brian

.. index::
   pair: example usage; NeuronGroup
   pair: example usage; run
   pair: example usage; SpikeGeneratorGroup
   pair: example usage; reinit
   pair: example usage; sqrt
   pair: example usage; Connection
   pair: example usage; SpikeCounter

.. _example-frompapers-computing with neural synchrony-coincidence detection and synchrony_Fig1F_ROC:

Example: Fig1F_ROC (frompapers/computing with neural synchrony/coincidence detection and synchrony)
===================================================================================================

Brette R (2012). Computing with neural synchrony. PLoS Comp Biol. 8(6): e1002561. doi:10.1371/journal.pcbi.1002561
------------------------------------------------------------------------------------------------------------------
Figure 1F.
(simulation takes about 10 mins)

Coincidence detection is seen as a signal detection problem, that of detecting a given depolarization in a background
of noise, within one characteristic time constant. The choice of the spike threshold implements a particular trade-off
between false alarms (depolarization was due to noise) and misses (depolarization was not detected).

Caption (Fig 1F). Receiver-operation characteristic (ROC) for one level of noise,
obtained by varying the threshold (black curve). The hit rate is the
probability that the neuron fires within one integration time constant t
when depolarized by Dv, and the false alarm rate is the firing
probability without depolarization. The corresponding theoretical
curve, with sensitivity index d' =Dv/sigma, is shown in red.

::

    from brian import *
    from scipy.special import erf
    
    def spike_probability(x): # firing probability for unit variance and zero mean, and threshold = x
        return .5*(1-erf(x/sqrt(2.)))
    
    tau_cd=5*ms     # membrane time constant (cd for coincidence detector)
    tau_n=tau_cd    # input is an Ornstein-Uhlenbeck process with the same time constant as the membrane time constant
    T=3*tau_n       # neurons are depolarized by w at regular intervales, T is the spacing
    Nspikes=10000   # number of input spikes
    T0=T*Nspikes    # initial period without inputs, to calculate the false alarm rate
    N=500           # number of neurons, each neuron has a different threshold between 0. and 3.
    w=1             # synaptic weight (depolarization)
    sigma=1.        # input noise s.d.
    sigmav=sigma*sqrt(tau_n/(tau_n+tau_cd)) # noise s.d. on the membrane potential
    print "d'=",1./sigmav # discriminability index
    
    # Integrate-and-fire neurons
    eqs='''
    dv/dt=(sigma*n-v)/tau_cd : 1
    dn/dt=-n/tau_n+(2/tau_n)**.5*xi : 1
    vt : 1 # spike threshold
    '''
    neurons=NeuronGroup(N,model=eqs,threshold='v>vt',reset='v=0',refractory=tau_cd)
    neurons.vt=linspace(0.,3,N) # spike threshold varies across neurons
    counter=SpikeCounter(neurons)
    
    # Inputs are regular spikes, starting at T0
    input=SpikeGeneratorGroup(1,[(0,n*T+T0) for n in range(Nspikes)])
    C=Connection(input,neurons,'v',weight=w)
    
    # Calculate the false alarm rate
    run(T0,report='text')
    FR=tau_cd*counter.count*1./T0
    # Calculate the hit rate
    counter.reinit()
    run(Nspikes*T,report='text')
    HR=counter.count*1./Nspikes-FR*(T-tau_cd)/tau_cd
    
    # Prediction based on Gaussian statistics
    FRpred=spike_probability(neurons.vt/sigmav)
    HRpred=spike_probability((neurons.vt-w)/sigmav)
    
    # Figure
    plot(FR*100,HR*100,'k')          # simulations
    plot(FRpred*100,HRpred*100,'r')  # theoretical predictions
    plot([0,100],[0,100],'k--')
    plot([0,100],[50,50],'k--')
    xlim(0,100)
    ylim(0,100)
    xlabel('False alarm rate (%)')
    ylabel('Hit rate (%)')
    
    show()